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Abstract 

A combined experimental and computational study aimed at high-resolution 3D 
imaging, visualization, and numerical reconstruction of fiber-reinforced polymer 
microstructures at the fiber length scale is presented. To this end, a sample of 
graphite/epoxy composite was imaged at sub-micron resolution using a 3D X-ray 
computed tomography microscope. Next, a novel segmentation algorithm was 
developed, based on concepts adopted from computer vision and multi-target 
tracking, to detect and estimate, with high accuracy, the position of individual fibers 
in a volume of the imaged composite. In the current implementation, the 
segmentation algorithm was based on Global Nearest Neighbor data-association 
architecture, a Kalman filter estimator, and several novel algorithms for virtual- 
fiber stitching, smoothing, and overlap removal. The segmentation algorithm was 
used on a sub-volume of the imaged composite, detecting 508 individual fibers. The 
segmentation data were qualitatively compared to the tomographic data, 
demonstrating high accuracy of the numerical reconstruction. Moreover, the data 
were used to quantify a) the relative distribution of individual-fiber cross sections 
within the imaged sub-volume, and b) the local fiber misorientation relative to the 
global fiber axis. Finally, the segmentation data were converted using commercially 
available finite element (FE) software to generate a detailed FE mesh of the 
composite volume. The methodology described herein demonstrates the feasibility 
of realizing an FE -based, virtual-testing framework for graphite/fiber composites at 
the constituent level. 

INTRODUCTION 

Three-dimensional imaging combined with numerical modeling of fiber 
reinforced polymer (FRP) composites at the constituent (i.e. fiber) length scale 
holds promise of significantly improving overall understanding of how these 
complex material systems deform, damage, and ultimately fail. This knowledge can 
eventually be incorporated into a multiscale “physical and virtual testing” 
framework [1-7], which will 1) enable discovery of new stochastic- and physics- 
based damage models, 2) reduce reliance on pure empiricism in design and 
certification of FRP structures [1,3], 3) enable development of revolutionary 
structural-life prognosis tools like the Airframe Digital Twin [8], and 4) enable 
accelerated discovery of new ultra-resilient multifunctional composites. 

In the past, 3D imaging and numerical modeling of “non-idealized” FRP 
microstructures has been restricted by reliance on surface-based imaging techniques 
(e.g. destructive sectioning coupled with optical and scanning-electron microscopy) 
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and by high computational costs associated with detailed 3D numerical simulations. 
However, rapid advances in non-destructive 3D micro- and nano-scale X-ray based 
imaging techniques 1 [9-14] and ongoing improvements in existing supercomputing 
hardware and numerical algorithms [8,15,16] are beginning to provide a means of 
overcoming such restrictions. 

Recently, a number of notable experimental studies have demonstrated the 
capability of high-resolution X-ray computed tomography (CT) techniques to 
visualize microstructure, micro-scale deformation, and damage in various fiber 
reinforced composites. For example, Latil et al. [17] used synchrotron radiation 
computed tomography (SRCT) to visualize in-situ 3D deformation of a 
fluorocarbon fiber bundle (with single fiber dia. of 150 pm) saturated with oil and 
subjected to transverse compression. In this study, image analysis was used to 
quantify position, orientation, displacement, and contact of individual fibers within 
the deforming bundle. In another study, Wright et al. [10] used sub-micron 
resolution SRCT to visualize post-mortem resin cracking and fiber rupture around a 
semi-circular notch in an aerospace-grade graphite/epoxy composite tensile coupon. 
In a later study, Wright et al. [12] utilized the same imaging technique to observe 
in-situ, micro-scale damage evolution in a double edge notched sample loaded in 
uniaxial tension. Perhaps the most impressive of these studies is the one conducted 
by Bale et al. [18], which utilized 1.3 pm resolution SRCT to visualize, in-situ, the 
sequence of microcracking in SiC-based composites subject to tensile loading at 
temperatures above 1600 °C. 

There are a number of other excellent examples of leveraging 2D and 3D 
imaging techniques (e.g. from X-ray CT) and numerical simulations across various 
length scales and material systems. For several examples of such approaches 
applied to composites, cellular materials, and biomaterials the reader is referred to a 
review article by Okereke et al. [7]. For recent examples of microstructural imaging 
of metallic materials using high-energy X-ray diffraction microscopy (HEDM) and 
supercomputer simulations based on HEDM data, the reader is referred to Refs. 
[19,20]. 

To the author’s knowledge, converting X-ray CT derived, constituent-level FRP 
data into high-fidelity, 3D numerical models has not yet been described in the 
literature. Until recently, most 2D and 3D microstructural finite element (FE) 
models have been created based on idealized or statistic ally-representative fiber 
arrangements (e.g. [21-30]), or based on 2D scanning-electron or optical 
microscopy image data (e.g. [25]). In the case of most 3D FE models, a 2D 
representation of the microstructure is simply extruded along the fiber-length 
direction to create the third dimension. While these approaches may prove very 
valuable in some cases (e.g. for modeling deformation and damage of relatively 
straight, unidirectional composites), they do not reflect the inherent variability in 
real FRP composites. Inclusion of 3D details, like precise fiber spacing, fiber 
waviness and misalignment, resin flaws, and foreign debris in numerical models 
may prove critical in fully realizing a multiscale virtual-testing framework [1]. 

Despite recent advances in high-resolution X-ray CT imaging techniques, 
conversion of X-ray CT-derived data into high-quality, realistic, numerical models 
of actual aerospace-grade FRP composites remains a significant challenge. One 
reason is that most available methods for converting 3D X-ray CT data (e.g. 3D 
stack of 2D gray-scale digital images) into numerical models (e.g. using gray-scale 
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thresholding) rely on a significant contrast in the attenuation coefficient between 
individual phases of the imaged material. However, for two-phase material systems 
like graphite/epoxy composites imaged at ~0.5- 1.0 pm resolution , the relatively 
low contrast between fibers and resin, small fiber diameter (-5-7 pm), varying 
cross-sectional size and shape of individual fibers (e.g. AS4, T300 fibers), and 
dense packing (fiber volume fraction -55-67%) render most phase-segmentation- 
based conversion codes inadequate. 

In light of the above, the work described herein and in the companion study [31] 
demonstrates a combined experimental and computational framework aimed at 
high-resolution 3D imaging, visualization, and numerical reconstruction of FRP 
microstructures at the fiber length-scale. To this end, this study adopts concepts 
from computer vision and multi-target tracking to detect and estimate, with high 
accuracy, the position of individual fibers in a volume of a graphite/epoxy 
composite imaged using a 3D X-ray microscope. 

In what follows, this manuscript begins with a description of the graphite/epoxy 
specimen followed by details of the 3D X-ray microscopy imaging. Next, the 
detection/estimation framework, including the template matching and fiber- 
tracking/assignment algorithms are presented, followed by a description of FE mesh 
generation. The manuscript is concluded by presenting results from the 
reconstruction, statistical quantification of the microstructure, and FE meshing. 

MATERIAL SYSTEM AND SPECIMEN PREPARATION 

The material system used in this study was a 24 ply, unidirectional, AS4/3501-6 
graphite/epoxy composite with a nominal single-ply thickness of 0.132 mm and 
average fiber diameter of 6.77 ± 0.28 pm. The specimen thickness was chosen to 
include at least 2-3 plies of the AS4/3501-6 composite while remaining thin enough 
to maintain quality of the scan and reduce scan time. The resulting “matchstick- 
type” specimen depicted in Fig. la had final dimensions of 0.46 mm x 0.33 nun x 
20.0 mm in the width (x), thickness (y), and length (z) directions, respectively. The 
z-direction corresponds to the fiber-length direction. Prior to X-ray CT scanning, 
the sample was adhered to a small finishing nail and mounted in the X-ray 
microscope as shown in Fig. lb. 

3D X-RAY MICROSCOPY 

The 3D X-ray microscopy was performed using an Xradia 500 Versa® 3 3D X- 
ray microscope at the Cornell University Biotechnology Resource Center Imaging 
Facility. The scan was perfonned over a 360° rotation using 4000 projections, 60 
kV voltage, 5 W power, 14 s exposure time, and 40X objective lens. The resulting 
voxel size was approximately 460 mn and the total scan time was 16 hours. 

Reconstruction of the attenuation data was performed using filtered back- 
projection, producing a stack of 949 cross-sectional, gray-scale digital images. The 
resulting cylindrical volume of the composite, depicted in Fig. 2, was 414 pm in 
diameter and 418 pm in height. For the purpose of subsequent algorithm and code 
development, a smaller, 169x169x169 pm sub-volume (approximate size of the 
cross-section shown with red dashed lines in Fig. 2) was extracted and used 
throughout this work. The sub-volume, the dimensions of which are approximately 
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50 times the average AS4-fiber radius, was chosen to be statistically representative 
of the bulk composite [23]. 


THE FIBER SEGMENTATION ALGORITHM 

The fiber segmentation algorithm described herein is divided into two steps, 
both implemented in Matlab® 4 . In the first step, a template-matching (TM) 
algorithm is used to detect (or “measure”) the 2D coordinates of individual-fiber 
centroids in each cross-sectional image that composes the volume shown in Fig. 2. 
In the second step, a multi-fiber assignment and tracking algorithm uses the TM 
detection data to detennine 3D coordinates of individual fibers within the imaged 
volume. A description of each step is presented in the following two sub-sections. 
For a more detailed discussion of the assignment and tracking algorithm, the reader 
is referred to a companion study by Whitacre and Czabaj [31]. 

Template Matching Algorithm 

For each 2D cross-sectional image in the 3D stack of tomographic images, the 
centroid coordinates of individual fibers are detected based on a TM algorithm [32, 
33]. This algorithm relies on computation of a normalized cross-correlation ( ncc ) 
score with the following expression: 


where (x,y) are the pixel coordinates within the cross-sectional image,/] is the MxN 
pixel cross-sectional image, fr is the mxn pixel template, f is the mean value of the 
mxn region within f) centered at (x.y), and f T is the mean value of the mxn template. 
The summation is perfonned over the mxn region for both fj and fr. Based on Eq. 
(1), a perfect correlation is obtained when the ncc score equals to 1, and perfect 
anti-correlation is obtained when the ncc score equals to -1. 

The process of detecting fiber centroids in a 2D cross-sectional image is 
described sequentially in Fig. 3. For clarity, a small, 70 x 70 pm subset of the entire 
cross-section (see Fig. 3a) is used in the current discussion. The TM procedure 
begins with the selection of several templates that are representative of the fiber 
cross-sections in the entire data set. Selection of individual templates was done by 
trial-and-error with the aim of maximizing the overall number of true detections, 
while minimizing the number of false ones. After several iterations, a total of six, 
17x17 pixel templates (shown in Fig. 3b), were selected and used throughout this 
study. Next, the ncc score was computed using Eqn. (1) for the entire image seen in 
Fig. 3a based on template #1. The resulting map of the ncc score is shown in Fig. 
3c. In this figure, the gray scale corresponds to ncc score between -1 (black) and 1 
(white). The map of the ncc score was then thresholded to values between 0.65-1.0 
to remove regions of low correlation. This results in a new image showing local 
peaks of ncc score, as seen in Fig. 3d. Next, a Matlab® function called 
imregionalmax was used to find centers of the local peaks in Fig. 3d. The resulting 
detections based on template #1 were superimposed onto the original image (Fig. 
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3a) using black crosses, and shown in Fig. 3e. Because a single template typically 
fails to detect all fiber cross-sections within any given image, the above process was 
repeated for the remaining five templates. In the final step, detections from all six 
templates were grouped together using a nearest-neighbor search algorithm with a 
search radius equal to one-half of the average fiber diameter. The resulting clusters 
of detections were then averaged, producing a final set of true and spurious 
detections, as seen in Fig. 3f. 

Fiber Assignment and Tracking Algorithm 

The TM algorithm generates a set of 2D detections across multiple images that 
can be used to accurately estimate the number and 3D coordinates of individual 
fibers throughout the imaged volume. Since the TM algorithm does not always 
detect all fiber cross-sections within each 2D image composing the 3D volume, and 
can sometimes generate spurious detections, a robust fiber estimation algorithm is 
required. 

To this end, a multi-fiber tracker based on a Kalman filter [34] is developed. 
The Kalman filter is a well established approach in target tracking, providing an 
optimal minimum mean square error estimate when the dynamics and 
measurements of the tracked object are linear and the errors are Gaussian. 

The Kalman filter is divided into predictor and corrector steps. In the prediction 
step, a fiber cross-section state (e.g. coordinates) in image k+1 is estimated based 
on the state from the previous image, k, using a linear extrapolation model. In the 
correction step, the new updated state estimate is generated based on an optimal 
correction using the corresponding TM detection at k+1. The assignment between 
the predicted state estimate and a corresponding TM detection (i.e. from among all 
possible TM detections at k+1 ) is achieved based on a Global Nearest Neighbor 
data association scheme [35]. 

In the current implementation of the fiber assignment and tracking algorithm, it 
is assumed that for each cross-sectional image the state of each fiber cross-section, 
x, can be described with a 2D position, P = [xy], and velocity, V= [x,y]. Here, V is 
assumed to be a function of spacing, r, between two consecutive tomographic 
images, equal to 460 nm (or 1 pixel). Further, it is assumed that all fibers are 
continuous and that no two fibers can overlap. 

The algorithm begins with the assignment of a unique fiber ID, i, to each 
detection within the first image, k = 0, in the stack. Next, for each detection i in 
image k, a new state estimate, x ik , is initialized: 

T 

%i,k ~ [xi,o ytfi 0 o] , (2) 

where jc/,o and yi t o are fiber centroid coordinates obtained from the TM algorithm at 
k =0. In Eqn. (2), it is assumed that each cross-section has zero initial velocity. 

After assignment of a unique ID and an initial state to each fiber at k, a new 
state estimate (or track), x ik+1 , and covariance matrix, P £;fc+1 , for the next image, 
k+1, are predicted using the standard Kalman filter prediction, 

Xi,k+i\k ~ ( 3 ) 

Pi,k+l\k = < &Pi,k\k < & T + Q- (4) 
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Note that the notation “a c \ b " represents an estimate of a in image c given 
observations up to, and including image b. A piecewise constant white acceleration 
model is used to simulate the change in fiber state between consecutive images 
[35]. This gives a state transition matrix, cD, as 


d> = 


h 

[ 0 2 



and the process noise covariance, Q, as 


Q — on 




(5) 

( 6 ) 


In Eqns. (5) and (6), h denotes a 2x2 identity matrix, O 2 is a 2x2 null matrix, and co 
is a process noise scaling factor. In Eqn. (3), the covariance matrix, Pi ik+1 \ k , is a 
measure of the estimated uncertainty of the new predicted state, x ik+1 \ k . 

After extrapolation, an assignment step makes the most likely assignment 
between each of the predicted tracks, x ik+1 , and the available detections za+i from 
TM data at k+1. In the current implementation, this assignment is based on the 
Munkres algorithm [36], requiring that it assigns, at most, one detection to a track, 
and no single detection to more than one track. Additional details of the assignment 
algorithm are given in [3 1]. 

Next, each track that is assigned to a new detection at k+1 is corrected with the 
standard Kalman fdter update, given by 


K ~ P i,k+l\kH T S i, k +\ 
%i,k+l\k+l — X iik+1 \k T Kv 
Pi,k+l\k+l Pi,k+l\k KH P t ; /c+i|/c? 


(V) 

( 8 ) 

(9) 


where K is the Kalman gain, H is the measurement matrix defined as H = [I 2 0 2 ], 
is the residual covariance defined as .S'/ k ■ / HP L k+1 \ k H r +R , R is the 

measurement error covariance (assumed constant for all measurements), and v is the 
measurement residual 


v = (Hx i]k+1{k -z iik+1 ). (10) 

The measurement residual measures the discrepancy between the predicted 
measurement, Hx ik+1 \ k and the assigned measurement Zt,k+i- 

For each track, the process of state prediction and correction described in Eqns. 
(3)-(10) is repeated over the entire stack of tomographic images. Throughout, any 
TM detections that did not get assigned to existing tracks are used to initiate new 
tracks. Additionally, tracks that have missed more than one update due to missed 
detection assignment are removed from further tracking [31]. The resulting set of 
tracks is next smoothed using a two-pass, forward and backward smoother [37]. 
Next, all tracks with length less than 169 pm (i.e. ones that do not span the entire 
volume) are compared to each other to determine which tracks correspond to the 
same fiber. Tracks that are determined to belong to the same fiber are then stitched 
together, and once more smoothed. The procedure of smoothing and stitching is 
repeated several times until there are no remaining fragmented tracks. 
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In the final step, all fibers are assigned an average fiber diameter of 6.77 pm 
and the entire data set is examined for possible overlap. In this study, fiber overlap 
results from errors associated with the TM algorithm, and from inherent variability 
in diameter of AS4 fibers. Any overlap between pairs of fibers is removed using 
two additional algorithms. The first algorithm finds fibers with significant overlap 
(i.e. > 20% of average fiber diameter) and keeps the most likely fibers without 
overlap. The second algorithm optimally adjusts positions and diameters of all 
overlapping fibers until all overlap is removed. A detailed description of smoothing, 
stitching, and overlap removal algorithms is further presented in [31]. 

FE MESH GENERATION 

After fiber segmentation, a multi-step algorithm was used to convert the 
resulting virtual-fiber data (i.e. coordinates, diameters) into a detailed 3D FE mesh 
of the composite. First, 3D coordinates of each virtual-fiber were imported into 
SolidWorks® 5 CAD software and converted into a 3D spline. Next, at the base of 
each spline, a circle with diameter corresponding to the respective fiber was 
sketched, and extruded along the spline to create a 3D solid model. The solid model 
generation was repeated for all segmented fibers using an automated Visual Basic 
macro routine. The resulting set of individual 3D solid-fibers was converted into a 
solid model assembly and subsequently imported into Abaqus/CAE® 5 using a 
binary parasolid file format. 

In Abaqus/CAE, a solid bounding box was created around the imported fibers, 
and a boundary-preserving subtraction/merge routine was used to create a solid 
region corresponding to the epoxy resin. The subtraction/merge routine enables 
generation of a coherent mesh between individual fibers and the surrounding resin. 
The resulting continuous volume, which includes fibers and resin, was seeded and 
meshed using 3D quadratic tetrahedral (C3D10) elements. 

RESULTS AND DISCUSSION 

The fiber segmentation algorithm described in the previous section was applied 
to a 169x169x169 pm sub-volume of the data shown in Fig. 2. The segmented sub- 
volume was composed from 367 gray-scale images with a z-direction spacing of 
460 mu. For this data set, the segmentation algorithm was run on a desktop 
workstation with a quad-core Intel Core i5-3470, 3.2 GHz CPU with 4 GB of RAM. 
The computational cost of the entire process including TM, fiber tracking, fiber 
smoothing, fiber stitching, and overlap removal was approximately 2 hours of wall- 
clock time. In addition to the automatic segmentation, a graphical-user interface 
(GUI) algorithm was developed and employed to manually segment 4 undetected 
fibers, and to remove 3 false tracks. 

Within the sub-volume of interest, the fiber segmentation algorithm identified a 
total of 508 individual fibers. Figure 4 shows a small, 70 x 70 pm, subset of the 
segmented volume at different z-coordinates with individual detections (numbered 
blue circles) overlaid on top of the tomographic data. As seen in Fig. 4, the position 
of most fiber cross-sections appears to be unchanged with the z-position, indicating 
good alignment with the z-axis. A few fibers (e.g. 302, 394, 509, 623) appear to be 
significantly misaligned from the z-axis, as indicated by the rapid change in the x-y 


5 SolidWorks and Abaqus/CAE are registered trademarks ofDassault Systemes. 
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coordinates of their respective centroids. Figure 5 shows the entire 169x169x169 
pm sub-volume of the imaged composite (Fig. 5a) alongside a solid model rendition 
of the corresponding fiber segmentation (Fig. 5b). Prior to visualization, the 
grayscale X-ray CT data in Fig. 5a were thresholded using Avizo 6 software to 
make the regions corresponding to the epoxy resin transparent. The segmented 
fibers in Fig. 5b were generated, assembled, and visualized in Solidworks® as 
described in the previous section. 

A qualitative (i.e. visual) examination of the data in Figs. 4 and 5 confirms that, 
within the resolution of the tomographic images, the fiber segmentation algorithm 
does a very good job of estimating shape and locations of individual fibers within 
the imaged volume. Additionally, the algorithm appears to accurately capture the 
size, shape, and distribution of the various resin rich regions. Currently, a more 
quantitative evaluation of the segmentation algorithm cannot be performed without 
higher resolution 3D data (e.g. from serial sectioning and optical imaging), which is 
currently a subject of ongoing research. 

After qualitative assessment of the segmentation algorithm, the 3D coordinates 
of the segmented fibers were used to quantify the relative distribution of individual 
fiber cross-sections within the imaged sub-volume. To this end, four x-y planes 
located at z = 0, 56, 113, and 169 pm were examined. The first measure of fiber 
distribution was based on computing the probability distribution function (PDF) of 
the nearest neighbor distance (NND). The NND is plotted in Fig. 6a as a function of 
spacing, h, between any two fibers, no/malized by the average AS4 fiber diameter, 
d ave . The second measure was based on computing a radial distribution function, 
g(h), which is related to the probability of finding a fiber cross-section within a 
circular annulus with internal radius, h, and thickness, Ah, centered around a 
reference fiber (shown in Fig. 6b). Several example applications of these two 
descriptors applied to 2D microstructural data of FRPs are presented in Refs. 
[26,28,38,39]. 

Examining both plots in Fig. 6, it is evident that most fibers within the imaged 
sub-volume are densely packed, with a majority of the fibers clustered within one 
average fiber-diameter of each other. Examining Fig. 6a, it appears that the 
probability of finding a nearest-neighbor fiber at h/d ave < 1 increases slightly with z, 
indicating greater fiber compaction at larger z. Qualitatively, this is consistent with 
the 3D data in Fig. 5, which shows a single diagonal fiber (#623) spanning the 
volume between z = 0 to 103 pm. In this region, the diagonal fiber effectively 
separates the volume into two distinct regions, decreasing the probability of NND 
for small z. Examining Fig. 6b, it appears that regardless of z, all fibers tend to 
cluster at h/d ave = 1, 2, 3, and 4. For h/d ave > 5, finding any two fibers away from the 
reference fiber is equally probable (i.e. g(h) — >1). 

The final quantitative examination of the microstructure was perfonned by 
computing an angle, a, which measures misorientation of any given fiber relative to 
the global z-axis. To simplify this calculation, each fiber was approximated as a 
straight line that intersects the fiber endpoints, and which is projected onto the x-z 
plane. The PDF of misorientation relative to the z-axis is presented in Fig. 7a as a 
function of a. For clarity, the a-axis in Fig. 7a was restricted to angles between ± 
5°. To help visualize the fiber z-axis misorientation, the same data were 
superimposed onto the 3D solid model of the segmented fibers using a color 


6 Avizo® is a registered trademark of Visualization Sciences Group, an FEI Company. 
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scheme as seen in Fig. 7b. In this figure, all fibers with a exceeding +5° were 
rendered with red color, while fibers with a exceeding -5° were rendered with dark 
blue color. As evident from Fig. 7, most of the fibers in the upper portion of the 
volume (i.e. above the red horizontal fiber #623) have slight negative misalignment, 
while most fibers in the lower portion have a slight positive misalignment. This 
difference in a between the two regions, along with the presence of a few 
significantly misaligned fibers in between, suggests that each sub-region 
corresponds to a distinct ply of the AS4/3501-6 laminate. 

In the final step, the segmented fibers depicted in Fig. 5b were imported into 
Abaqus/CAE® and used to create an FE mesh of the entire composite using the 
procedure outlined in the previous section. To generate a well-structured, high- 
quality mesh, a seed density of about 12 elements per fiber perimeter was used 
throughout. When necessary, Abaqus/CAE® was allowed to adaptively relax the 
specified element size constraint to account for complex geometrical features (e.g. 
close spacing and contact between individual fibers). Meshing was performed on a 
desktop workstation with a quad-core Intel Core i7-3720QM, 3.6 GHz CPU with 32 
GB of RAM, with a computation cost of approximately 1 wall-clock hour. The final 
FE mesh shown in Fig. 8, contained ~6.7 million quadratic, tetrahedral C3D10 
elements and ~21 million degrees of freedom (DOF). Considering the scope of this 
study, no attempts were made to perform a detailed mesh convergence study. 
However, to explore tractability of the current model, a static, linear-elastic 
simulation of uniaxial tension (along y-axis) was performed with FEAWD 7 parallel 
FE solver using 1000 Intel Ivy Bridge cores of NASA’s Pleiades supercomputer 
(total 184,800 cores available). The pre- and post-processing of the FE model and 
the results took several hours; however, the actual static simulation was completed 
in approximately 140 wall-clock seconds. The same simulation was completed in 
55 seconds using 5000 cores. 

The combined imaging, segmentation, and meshing procedure described herein 
provides a significant step towards realizing the idea of virtual testing of FRPs at 
the constituent scale. However, despite the very promising results shown herein, a 
number of important technical challenges remain. For instance, in the present study 
a relatively small, unidirectional matchstick specimen was used to obtain high 
quality X-ray CT data. For practical applications it may be necessary to acquire 
such data from larger, macro-scale components with off-axis fiber stacking 
sequences, over limited scan angles (i.e. <180°). As a result, the present 
segmentation algorithm must be extended to enable processing larger volumes of 
noisier and perhaps incomplete data. Finally, despite the relatively low computation 
cost (for a Pleiades-type supercomputer) of performing static, linear-elastic 
simulations, future inclusion of non-linear effects and damage (e.g. plasticity, 
fracture) into these models may increase the computational cost dramatically. With 
this in mind, new, adaptive-meshing algorithms that generate efficient meshes will 
have to be developed. 

CONCLUSION 

This manuscript describes a combined experimental and computational study 
aimed at high-resolution 3D imaging, visualization, and numerical reconstruction of 
fiber reinforced polymer microstructures at the fiber length-scale. To this end, a 


7 FEAWD was developed at the Cornell University Fracture Group (www.cfg.cornell.edu). 
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small sample of graphite/epoxy composite was imaged at sub-micron resolution 
using a 3D X-ray computed tomography microscope. Next, a novel segmentation 
algorithm was developed based on concepts adopted from computer vision and 
multi-target tracking to detect and estimate, with high accuracy, the position of 
individual fibers in a volume of the imaged composite. In the current 
implementation, the segmentation algorithm was based on a Global Nearest 
Neighbor data-association architecture, a Kalman filter estimator, and several novel 
algorithms for fiber stitching, smoothing, and overlap removal. The segmentation 
algorithm was demonstrated on a 169x169x169 pm sub-volume of the imaged 
composite, detecting 508 individual fibers. The segmentation data were 
qualitatively compared to the tomographic data, demonstrating high accuracy of the 
numerical reconstruction. Moreover, the data were used to quantify a) the relative 
distribution of individual fibers cross-sections within the imaged sub-volume, and 
b) the local fiber misorientation relative to the global fiber axis. Finally, the 
segmentation data were converted using commercially available finite element (FE) 
software to generate a detailed finite element mesh of the composite volume. The 
methodology described herein demonstrates the feasibility of realizing an FE-based 
virtual-testing framework for graphite/fiber composites at the constituent level. 
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Fig. 1 The AS4/3501-6 specimen a) prior to scanning and b) mounted inside Xradia 

500 Versa® microscope. 
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Fig. 2 AS4/3501-6 graphite/epoxy specimen: a) reconstructed volume and b) cross- 
sectional image with sub-volume selected for numerical reconstruction outlined in 

red dashed box. 




Fig. 3 Sequence of images describing the TM algorithm: a) 2D X-ray CT image, b) 
six fiber templates, c) map of ncc score from template #1, d) map of ncc score after 
thresholding, e) detections from template #1, f) detections from templates #1-6. 
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(a) 




(c) 


Fig. 4 Example of tracked fibers (blue circles) overlaid on top of the X-ray CT data 
at a) z = 0 pm, b) z = 56 pm, c) z =1 13 pm, and d) z = 169 pm. 




(a) (b) 

Fig. 5 3D rendition of (a) the X-ray CT data, and (b) solid model representation of 

the segmented fibers. 
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ig. 6 Statistical description of fiber arrangement at different z-coordinates based 
on a) the nearest neighbor distance and b) radial distribution function. 



Fig. 7 (a) PDF of fiber misalignment relative to the z-axis, and (b) 3D CAD 
rendition of the segmented fibers with color scheme representing the misalignment 

angle, a. 
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Fig. 8 Finite element mesh created using A baq us/CAE (-6.7M tetrahedral 

elements, ~27M DOF). 
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